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ABSTRACT 

We used the N-body code of iHernquist &: Ostrikerl (|1992[ ) to build a dozen cuspy 
(7 ~ 1) triaxial models of stellar systems through dissipationless collapses of initially 
spherical distributions of 10 6 particles. We chose four sets of initial conditions that re- 
sulted in models morphologically resembling E2, E3, E4 and E5 galaxies, respectively. 
Within each set, three different seed numbers were selected for the random number 
generator used to create the initial conditions, so that the three models of each set 
are statistically equivalent. We checked the stability of our models using the values of 
their central densities and of their moments of inertia, which turned out to be very 
constant indeed. The changes of those values were all less than 3 per cent over one 
Hubble time and, moreover, we show that the most likely cause of those changes are 
relaxation effects in the numerical code. We computed the six Lyapunov exponents 
of nearly 5,000 orbits in each model in order to recognize regular, partially and fully 
chaotic orbits. All the models turned out to be highly chaotic, with less than 25 per 
cent of their orbits being regular. We conclude that it is quite possible to obtain cuspy 
triaxial stellar models that contain large fractions of chaotic orbits and are highly 
stable. The difficulty to build such models with the method of ISchwarzschildl (|1979T) 
should be attributed to the method itself and not to physical causes. 

Key words: Galaxies: elliptical and lenticular, cD - Galaxies: kinematics and dy- 
namics - methods: numerical - Physical data and processes: chaos. 



1 INTRODUCTION 

Self-consistent models of spherical and disklike stellar sys- 
tems are relatively simple to build u sing standard textbook 
methods (jBinnev fc Tremainell2008l ). but special techniques 
are necessary to obtain triaxial models adequate to repre- 
sent elliptical galaxie s. One of the most p opular of those 
techniques is due to ISchwarzschildl (1 19791 ): one chooses a 
reasonable mass distribution for the system, obtains the 
corresponding potential and computes a large library of or- 
bits in that potential selecting suitable initial conditions; 
weights are then assigned, according to the time spent on 
each orbit and in every region of space, and used to set 
up a system of linear equations linking the mass density in 
each region to the fractions of the different types of orbits, 
which can be finally obtained solving the system. A differ- 
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ent method, which we will refer to as the N-body method, 
was due to ISparke fc Sellwoodl l|l987h : adopting an initial 
distribution of mass points, one integrates their equations 
of motion with an N-body code until an equilibrium con- 
figuration is reached; the potential is then fixed and fit- 
ted with an adequate smooth approximation that allows 
the computation of orbits using the positions and veloci- 
ties of the bodies as initial conditions. The initial distri- 
bution of bodies can be chosen in different ways to obtain 
a triaxial system: slow deforma tion of a spherical system 
(|Hollev-Bockelmann et al.ll200l | ). cold collapse of a spher- 
ical system dVoglis, Kalapotharakos fc Stavropoulosl |2002| ; 
M uzzio. Carpintero fc Wachlinll2005l) . merging of stellar sys- 
tems ( Jesseit. Naab fc Burkertll2005f ). and so on. Neverthe- 
less, with the N-body method it is not possible to have 
the degree of control over the final configuration that one 
has with Schwarzschild's method. In the end, both meth- 
ods yield the same result (i.e., a self-consistent stellar sys- 
tem and an analysis of its orbital content), but the method 
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of lSchwarzschndl (| 19791 ) starts with a chosen mass distribu- 
tion, obtains orbits and finds which fractions of those orbits 
are needed to have a self-consistent system with that mass 
distribution, while the N-body method begins obtaining a 
self-consistent system and then performs the orbital analy- 
sis. 

Chaotic moti ons are frequent i n st ellar systems 
(|ContopouIosll2004D , and lSchwarzschildl (| 19931 ) realized that 
triaxial systems should include chaotic orbits but, when 
he introduced them in his models, the models evolved 
on time scales of the order of a Hubble time and were 
not truly stable. This problem was aggravated after it 
became clear that cuspy models (where near the cen- 
tre the density, p(r), is proportional to r where r is 
the radi us and 1 ^ 7 ^2) w ere needed for elliptical 
galaxies (jMerritt fc Fridmanlll996l ). Nevertheless, it is per- 
fectly possible to obtain stable triaxial models, even cuspy 
ones, that contain l arge fractions of chaotic orbits using 
the N-body method dVoglis et alj|2002l; iMuzzio et al.ll200a 
lAquilano et al.ll2007l ; lMuzzio, Navone fc Zorzill2009l ). so that 
several of those authors argued t hat the difficulty to build 
such systems with the method of ISchwarzschildl (| 19T9h was 
due to that method itself and not to a physical cause. 

Clearly, chaot ic orbits are dif f icult t o deal with when us- 
ing the method of ISchwarzschildl l|l993h because they cover 
different regions at different times. Worse still, in cuspy tri- 
axial potentials chaotic orbits tend to be extremely sticky 
l|Siopis fc Kandrupll2000l ) . and sticky orbits behave like reg- 
ular orbits for long periods of time and chaotically at other 
intervals. Let us assume, as an example, that chaotic orbits 
occupy part of the time an elongated region of space and an- 
other part of the time a nearly spherical region. Therefore, 
the orbits that occupy mostly an elongated region during 
the integration time used to prepare the library of orbits 
will have large weights in that region and low weights out- 
side it; conversely, the orbits that occupy mostly the nearly 
spherical region during the same integration time will have 
a more or less even distribution of weight values over this 
region. Now, when these weights are used to establish the 
fractions of orbits that make up the triaxial system, the or- 
bits of the former group will be strongly fav oured and the 
equilib rium system yielded by the method of ISchwarzschildl 
i| 19791) will include many of them and few of the other group. 
If we then let the system evolve, the former orbits will tend 
to fill in a more spherical region as time goes by and, con- 
versely, the latter orbits will adopt a flatter distribution. 
But, since the model included less orbits that originally 
had a more spherical distribution, t he model will become 
rounder, as shown by the Table 6 of ISchwarzschildl (|l993l ) 
fo r most of his models. The use o f longer integration times, 
as lCapuzzo-Dolcetta et all (|2007l ) advocated, does not guar- 
antee the success of the method either, because the average 
behaviour of a chaotic orbit over a certain interval, no matter 
how long, does not necessarily coincide with its behaviour 
over a different interval, or even a subinterval, of the inte- 
gration time. To make matters worse, the number of orbits 
typically computed for Schwarzschild's method is not very 
large, from several hundreds in the works of Schwarzschild 
himself up to about 10, 000 in more recent work like that of 
ICapuzzo-Dolcetta et al.l (|2007l ) , so that they do not provide 
a strong enough statistical basis. 

On the other hand, the N-body method uses large num- 



bers of bodies (of the order of 10 6 in our own more recent 
works) guaranteeing good statistics and the model is created 
by the evolution of the system itself, that is, the dynamics 
take care of favoring certain orbits at the expense of others 
so as to reach a self-consistent state. Is that state a sta- 
ble one? If all the orbits are regular it certainly is because, 
once reached that state, the orbits will continue filling in the 
same regions of space for ever, in other words, we will have a 
static equilibrium. If the model includes chaotic orbits, how- 
ever, it is conceivable that, once in equilibrium, when part 
of the chaotic orbits evolves to occupy a rounder space re- 
gion, another part of the chaotic orbits which had a rounder 
distribution evolves towards a more elongated one and fill 
in the phase space regions left vacant by the former, that 
is, we might have a dynamic rather than static equilibrium. 
While all this is concievable there is no guarantee that it will 
actually happen, but the past investigations that resulted in 
stable models that contained high fractions of chaotic orbits 
undoubtedly support this view. 

Therefore, we want to investigate whether it is possible 
to obtain self-consistent models of cuspy triaxial sytems that 
are stable. That is, we will deal with the type of model that 
is most likely to contain high fractions of chaotic orbits and 
most difficult to obtain with Schwarzschild's method. We 
have already obtain ed such models with the N-body method 
(|Muzzio et alj|2009l ) but, in order to compensate for t he soft- 
ening needed by the N-body code of L.A. Aguilar (| White! 
ll983l : lAguilar fc Merrittlll990l ) we had to introduce an addi- 
tional potential. Here we wi ll use the code of L. Hernquist 
|Hernquist fc Qstrikerl Il992l ) that does not need softening 
and that uses an expansi on of the potential based on the 
model of lHernquistl (|l990h that is particularly adequate for 
cuspy models. 

The next section gives a description of how we obtained 
our models and their main properties. Section 3 presents 
our results on the stability and chaoticity of the models and 
Section 4 summarizes our conclusions. 



2 MODELS OF CUSPY TRIAXIAL STELLAR 
SYSTEMS 

2.1 Model building 

We built our mod els following the recipe of 
lAguilar fc Merrittl (|l990lV_ just as we have done in our pre- 
vious i nvestigations (see Muzzio et al.ll2005l ; lAquilano et al.l 
120071 ; M 

uzzio et al.l |2009| , for example): we randomly 
created a spherical distribution of 10 6 particles with a 
density distribution inversely proportional to the distance 
to the centre and a Gaussian velocity distribution, and 
we let it collapse follow i ng th e evolution with the code 
of iHernouist fc Qstrikerl (|l992T l; due to the radial orbit 
instability, the result is a triaxial system. The gravitational 
constant, G, the radius of the sphere and the total mass 
are all set equal to 1 and the collapse time (the one needed 
to reach the maximum potential energy) turns out to be 
also very close to unity. We followed the initial collapse for 
7 time units (t.u., hereafter) and, then, we eliminated the 
particles with positive energy, we determined the principal 
axes of the inertia tensor of the 80 per cent most tightly 
bound particles and we rotated the system so that its major, 
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intermediate and minor axes coincide, respectively, with 
the x, y and z axes of coordinates. Subsequently, we let the 
system evolve for another 300 crossing times (T cr , hereafter, 
see Table 1 below), eliminated the less bound particles that 
had not yet reached equilibrium (less than 2 or 3 per cent 
of the particles in all cases) and aligned it again with the 
system of coordinates. We then performed a last run of 600 
T cr to, first, obtain the final model after the initial 100 T cr 
and, second, check its stability over the final 500 T cr . We 
must recall that we are only interested in obtaining models 
morphologically similar to elliptical galaxies and not in 
creating them through a realistic process. Besides, we want 
to make sure that our models have reached the equilibrium 
state before testing them for possible evolutionary effects. 
As we will show later, one Hubble time is of the order of 200 
T cr for our models, so that the integration times used here 
are indeed much larger than that time and, as indicated, do 
not correspond to a realistic process of galaxy formation. 

After several trials, we chose an integration step of 
0.0025T cr which yielded integration errors in the energies 
of the individual bodies negligibly small compared to the 
chan ges in those energies due to the relaxation effects of the 
code (|Hernquist fe Barnesj|l990l ). With that choice, the total 
energy is conserved within about 0.1 per cent during most 
of the evolutions, except during the initial collapse when it 
is conserved only within about 1 or 2 per cent. 

iKalapotharakos. Efthvmiopoulos fc Vogiisl (|2008f ) in- 
vestig ated the approx imation of N-body realizations of mod- 
els of iDehnenl (|l993l) for ^ 7 1 with a g eneraliza- 
tion of the method of lHernquist fc Ostrikerl (| 19921 ) and they 
found that the choice of the radial basis functions seri- 
ously affects the results of the fractions of chaotic orbits 
or the distribution of the Lya p unov characteristic expo- 
nents. The model of iHernouistl Jl990l) corres ponds to the 
7 = 1 cas e of the models of Dchncnl (| 19930 so that the 
method of lHernquist fc Ostrikerl (1 19921 ). based on the for- 
mer, should be expected to yield good results for models 
wi th 7 = 1 and, in fact, t hat is one conclusion of the work 
of IKalapotharakos et al.l |2008l ). Therefore, it is important 
that we stick to models with 7 = 1, not only because they 
are fairly cuspy, but also because the N-body method we 
are using works best for such models. Notice that even our 
initial distribution (an sphere with density proportional to 
the inverse of the radius) obeys that 7 value. 

The code of iHernouist fc Ostrikerl (| 19921 ) allows one 
to choose the number of terms in the angular and radial 
expansions (l ma x and n max , respectively), so that we per- 
formed tests with different numbers of terms in the range s 
2 l max ^ 6 and 4 n max s: 10 |Zorzi fc Muzzidl2009D . 
For every (I 

max 1 Tlmax 

) pair we run three simulations that 
differed only in the seed number used to randomly gener- 
ate the initial positions and velocities, i.e., the three were 
statistically equivalent. The pairs (3, 6), (3.8), (3, 10), (4, 7), 
(4, 10), (5, 8) and (6, 8) yielded models that were only moder- 
ately cuspy, with 7 ~ 0.6 — 0.9, so that they were discarded. 
It might seem odd that, despite the use of a basis of ra- 
dial functions derived from the Hernquist model, which has 
7=1, the obtention of collapse models with such cuspiness 
is not guaranteed but, although the zero order radial func- 
tion is cuspy with 7 = 1, higher order terms are not and in all 
likelihood their contribution flattens the cusp. Most of the 
remaining (l ma x, n ma x) pairs might have been acceptable, as 



the reasons to prefer one to another were not as compelling 
as the low 7 values that made us reject those mentioned 
before. Pairs (2, 8) and (4, 5) yielded 7 values that departed 
moderately (less than 0.1) from 1, and the results from their 
different statistical realizations exhibited larger dispersions 
than other pairs. Of the three last pairs, all with lmax = 4, 
the one with n ma x = 6 showed the most consistent results 
among the different statistical realizations, so that we finally 
decided to use those numbers of terms in our expansions. 

2.2 The models 

Following lAguilar fc MerrittJ (|l990l ) we took the square roots 
of the mean square values of coordinates of the 80 per cent 
most tightly bound particles as the semiaxes of the sys- 
tem; the major, intermediate and minor axes will be dubbed 
a, b and c, respectively, hereafter. We built several models 
adopting different values of the dispersion for the velocity 
distribution, obtaining less elongated models for larger dis- 
persion values. Finally, we adopted four dispersion values 
that resulted in models with c/a values close to 0.8, 0.7, 0.6 
and 0.5, respectively (i.e., corresponding to elliptical galax- 
ies between Hubble types E2 and E5). The velocity disper- 
sion for the E5 model was essentially zero, so that more 
elongated models could not be obtained. A similar result 
had been obtained in our previous work (|Aouilano et al.l 
120071) . where we could not obtain models more flattened 
than E6 and we suggested that this might hint that merg- 
ers, rather than collapses, are needed to obtain the most 
elongated ellipticals. For each selected value of the veloc- 
ity dispersion, three different models (dubbed a, b and c) 
were obtained using three different seed numbers for the 
random number generator when creating the initial distri- 
bution of particles. In brief, we prepared a grand total of 
12 models, divided in 4 groups; the models in each group 
resemble E2, E3, E4 and E5 galaxies, respectively, and the 
three models within each group differ from each other only 
from a microscopic point of view, being essentially identical 
from the point of view of their macroscopic properties. As 
m our previo us work (|Muzziol 120061 ; lAquilano et~al1 120071 ; 
M UZZIO et al.ll2009l ), several models display figure rotation 
around the minor axis, i.e., they rotate very slowly even 
though their total angular momentum is zero. 

Table 1 summarizes the global properties of our mod- 
els: mass, crossing time, effective radius, central radial ve- 
locity dispersion, triaxiality, 7 and angular velocity of figure 
rotation. The effective radius was obtained from the (x, z) 
projection and, accordingly, the central radial velocity dis- 
persion was computed from the y components of the ve- 
locities of the 10, 000 particles closer to the centre on that 
projection. Triaxiality was evaluated from the semiaxes ob- 
tained from the 80 per cent most tightly bound particles 
as T = (a 2 — b 2 )/(a 2 — c 2 ). 7 was obtained as the slope of 
the logp(r) vs. log(r) line for the innermost 10,000 parti- 
cles binned in 100 particle bins. The angular velocity was 
computed from the angles formed by the major axis with its 
original position at different times of the final 500 T cr in- 
tegration. The Table clearly shows that, for a given model, 
the different realizations obtained changing the seed number 
in the random number generator have essentially the same 
global properties. The only exception is the angular veloc- 
ity which displays significant differences, particularly among 



4 A. F. Zorzi and J. C. Muzzio 



Table 1. Properties of our models. 



Model 


Mass 


T 




(7() 


T 




to 


E2a 
E2b 
E2c 


0.990 
0.990 
0.990 


0.721 
0.723 
0.721 


0.250 
0.253 
0.250 


0.894 
0.897 
0.902 


0.73 
0.78 
0.70 


1.052 ± 0.021 
0.997 ± 0.021 
1.006 ± 0.021 


0.0003 ± 0.0001 
0.0012 ± 0.0004 
0.0016 ± 0.0004 


E3a 
E3b 
E3c 


0.979 
0.979 
0.978 


0.659 
0.659 
0.657 


0.220 
0.221 
0.219 


0.934 
0.933 
0.929 


0.65 
0.68 
0.66 


1.056 ± 0.021 
1.062 ± 0.021 
0.987 ± 0.021 


-0.0004 ± 0.0001 
0.0002 ± 0.0001 
0.0010 ± 0.0001 


E4a 
E4b 
E4c 


0.911 
0.905 
0.908 


0.526 
0.518 
0.521 


0.182 
0.180 
0.180 


0.957 
0.974 
0.963 


0.62 
0.63 
0.62 


1.053 ± 0.019 
1.027 ± 0.021 
1.073 ± 0.022 


0.0001 ± 0.0003 
-0.0002 ± 0.0001 
-0.0002 ± 0.0002 


E5a 
E5b 
E5c 


0.906 
0.908 
0.907 


0.463 
0.466 
0.465 


0.159 
0.160 
0.158 


0.997 
0.986 
0.970 


0.46 
0.48 
0.46 


0.985 ± 0.022 
1.000 ± 0.019 
1.017 ± 0.020 


0.0020 ± 0.0007 
0.0059 ± 0.0004 
0.0164 ± 0.0003 



models E5a, b and c. Nevertheless, rotation is very slow in 
all cases and, except for the E5c case, lower than the 0.00975 
value of the investigation bv lMuzziq l|2006h . which revealed 
only very small differences between the fractions of chaotic 
orbits in the rotating and non-rotating system. The present 
E4 and E5 models are much more triaxial than the corre- 
sponding non-cuspy models of lAquilano et all (|2007t ) (0.98 
and 0.81, respectively), and the present mod el E4 is much 
more triaxial than the cuspy model E4c of M uzzio et al.l 
<|2009h (0.91). 

As in ou r prev ious work (see lAquilano et al.l 120071 ; 
iMuzzio et al.l 120091 . for exam ple), we chose galax- 
ies NGC1379 and N GC4697 (|Napolitano et all 120051 : 
iForbes fc Ponmanl IT999T ). whose mass-to-light ratio gradi- 
ents are zero, to obtain some estimate of the equivalence 
between our units and those of real galaxies, Comparing 
their observed values of R e (2.5 and 5.7kpc, respectively) 
and cto (128 and 180 fcms -1 , respectively) with those from 
our Table [ll we conclude that values between about 10 and 
36fcpc can be used as our length unit and values between 
about 0.07 and 0.20Gy as our time unit. Then, the Hubble 
time can be estimated as between 66 and 190 t.u., and we 
will adopt a value of 100 t.u., hereafter. 

We obtained the projected d istributions of the particle s 
on the (xz) plane and, following lAguilar fc Merrittl (|l990h . 
we adopted fixed axial ratios equal to those of the 80 per 
cent most tightly bound particles for each model (see be- 
low), and computed the surface density of particles in el- 
liptical shells containing 10, 000 particles per bin. Just as 
lAguilar fc Merrittl (|l99C ) had found for their models, ours 



follow very closely a de Vaucouleurs law and, again, differ- 
ences among the different random realization of the same 
model are negligibly small. As an example, Figure Q] shows 
the results for our models E5. 

Figure[2]shows, as an example, the logarithm of the den- 
sity versus the logarithm of the radius for the central part of 
the E5 models; different symbols correspond to models with 
different seed numbers and a straight line with slope 7 = 1 
is also shown for comparison. In this case the density was 
computed in spherical shells containing 100 particles each, 
so that the Poisson relative error of each point on the graph 
is about 0.043. Besides, Figure [3] presents, for the innermost 
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Figure 1. Surface density versus ij 1 / 4 plot for the E5 models 



region of model E5c, the logarithm of the density versus the 
logarithm of the radius for both the adopted model (circles) 
as for that same model evolved 100 t.u., or about a Hubble 
time, (plus signs) and we note that the slope of the central 
cusp is very well conserved. Similar results were obtained for 
the other models. 

Table [2] gives for each model the values of the major 
semiaxis and of the axial ratios for the 20, 40, ...and 100 
per cent most tightly bound particles. Again, the results 
obtained changing the randomly generated inital conditions 
for a given model are essentially the same. As in our previous 
models, there is a general trend towards larger axial ratios 
at larger distances from the centre, but the trend is weaker 
for the more flatened models and breaks closer to the centre, 
probably due to the influence of the cusp. 
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Table 2. Major semiaxes and axial ratios of our models 



Model 


Pvopevty 


20% 


40% 


60% 


80% 


100% 


E2a 


a 
b/a 
c/ a 


0.068 
0.754 
0.596 


0.117 
0.789 
0.684 


0.172 
0.837 
0.768 


0.300 
0.877 
0.826 


0.578 
0.930 
0.909 


E2b 


a 
b/a 
c/a 


0.068 
0.742 
0.602 


0.118 
0.777 
0.687 


0.175 
0.820 
0.763 


0.297 
0.870 
0.829 


0.574 
0.925 
0.906 


E2c 


a 
b/a 


0.067 
0.767 
0.605 


0.116 
0.801 
0.691 


0.173 
0.845 
0.769 


0.294 
0.882 
0.826 


0.578 
0.937 
0.913 


E3a 


a 
b/a 
c/a 


0.062 
0.753 
0.597 


0.115 
0.723 
0.564 


0.164 
0.781 
0.642 


0.294 
0.814 
0.694 


1.097 
0.852 
0.810 


E3b 


a 
b/a 
c/a 


0.065 
0.710 
0.554 


0.117 
0.701 

0.555 


0.165 
0.761 
0.634 


0.294 
0.802 
0.692 


1.139 
0.840 
0.803 


E3c 


a 
b/a 


0.065 
0.714 
0.559 


0.115 
0.715 
0.559 


0.165 
0.776 
0.635 


0.294 
0.808 
0.689 


1.144 
0.852 
0.807 


E4a 


a 
b/a 
c/a 


0.056 
0.731 
0.581 


0.104 
0.697 
0.516 


0.152 
0.735 
0.546 


0.233 
0.769 
0.581 


1.146 
0.775 
0.674 


E4b 


a 
b/a 
c/ a 


0.056 
0.728 
0.575 


0.102 
0.698 
0.517 


0.150 
0.733 
0.549 


0.231 
0.764 
0.583 


1.160 
0.792 
0.706 


E4c 


a 
b/a 
c/a 


0.055 
0.749 
0.588 


0.103 
0.706 
0.515 


0.152 
0.737 
0.539 


0.232 
0.765 
0.575 


1.180 
0.797 
0.708 


E5a 


a 
b/a 
c/a 


0.052 
0.824 
0.557 


0.095 
0.809 
0.508 


0.148 
0.815 
0.504 


0.216 
0.814 
0.515 


0.482 
0.893 
0.563 


E5b 


a 
b/a 
c/ a 


0.051 
0.846 
0.596 


0.095 
0.810 
0.509 


0.147 
0.802 
0.506 


0.222 
0.803 
0.508 


0.479 
0.879 
0.555 


E5c 


a 
b/a 
c/a 


0.051 
0.847 
0.583 


0.094 
0.825 
0.512 


0.145 
0.810 
0.500 


0.222 
0.810 
0.506 


0.496 
0.883 
0.569 



3 RESULTS AND ANALYSIS 
3.1 Stability 

We used the results of the 500 T cr long final evolution to 
check the stability of the central density and the semiaxes 
of the stellar systems. For this purpose, we computed at 
100 T cr intervals the central density from the 10, 000 parti- 
cles closer to the centre of each system, and the moments 
of inertia of the 80 per cent most tightly bound particles. 
These quantities changed almost linearly with time, so that 
we obtained their variations over a Hubble time from the 
corresponding best fitting straight lines and the results are 
presented in Table [3] 

Although most of the variations of the central density 




log(radius) 



Figure 2. Density versus radius plot for the E5 models 



4 | , , , , 1 , , , , 1 , , , , 1 , , , r 



o Original 

+ 1 Hubble time 




Figure 3. Density versus radius plot for the E5c model at differ- 
ent times. 

are not significant at the 3<r level it is suggestive that, except 
for that of model E5b (also not significant), they are all 
negative. Almost all the variations of the moments of inertia 
are highly significant and indicate a general decrease of the 
major axes as well as general increases of the intermediate 
and minor axes. But the changes are very small indeed, all 
of them being smaller than 3 per cent over one Hubble time. 

Interestingly, both the amounts and the senses of the 
changes are similar to those found i n our previous works 
(|Aquilano et al.1 120071 ; iMuzzio et"aH 120091 ). where we at- 
tribut ed them mainly to relaxati on effects of the multipolar 
code (|Hernquist fc Barnes! 1 1990h . Since in those works we 
had used the N-body code of Aguilar and we are now using 
that of Hernquist and Ostriker, we decided to perform the 
same checks we had made before using the models E5, i.e., 
those that exhibit the largest changes. The first check was to 
eliminate self-consistency, letting the systems evolve again 
for 500 T cr but keeping the coefficients of the expansion of 
the potential fixed at their initial values, and the results 
are shown in Table [4] The central density changes, although 
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Table 3. Percentage variations over one Hubble time. 



Model Dens.(%) Xmom.in.{%) Ymom.in.(%) Zmom.in.(%) 



E2a 


-0.61 


± 


0.29 


-0 


.63 


± 


0.09 





.64 


± 


0.13 





.20 


± 





09 


E2b 


-0.36 


± 


0.43 


-0 


.18 


± 


0.07 





.51 


± 


0.08 





■11 


± 





11 


E2c 


-0.83 


± 


0.59 


-0 


.45 


± 


0.10 





19 


± 


0.10 





.40 


± 





13 


E3a 


-0.76 


± 


0.53 


-1 


01 


± 


0.05 


1 


.36 


± 


0.11 





.78 


± 





10 


E3b 


-0.62 


± 


0.15 


-1 


.08 


± 


0.07 


1 


11 


± 


0.09 





.86 


± 


0. 


15 


E3c 


-0.56 


± 


0.18 


-1 


.04 


± 


0.11 


1 


.57 


± 


0.16 





.87 


± 


0. 


11 


E4a 


-1.05 


± 


0.65 


-1 


12 


± 


0.06 


1 


.16 


± 


0.04 


1 


.68 


± 





12 


E4b 


-0.86 


± 


0.86 


-1 


.21 


± 


0.09 


1 


19 


± 


0.08 


1 


.65 


± 





17 


E4c 


-0.73 


± 


0.56 


-1 


20 


± 


0.10 


1 


.36 


± 


0.23 


1 


70 


± 





18 


E5a 


-1.07 


± 


0.78 


-2 


17 


± 


0.11 


2 


12 


± 


0.14 


2 


.65 


± 





21 


E5b 


0.92 ± 


0.83 


-1 


96 


± 


0.13 


2 


00 


± 


0.08 


1 


.95 


± 





08 


E5c 


-1.23 ± 0.17 


-2 


02 


± 


0.09 


1 


73 


± 


0.21 


2 


.31 


± 


0. 


06 



somewhat lower than the corresponding values of Table [3l 
are again not significant at the 3cr level, but the changes of 
the moments of inertia are much lower and generally sig- 
nificant at the same level. This is what could be expected 
from changes due to relaxation effects of the N-body code, 
because they would be suppressed when turning off the self- 
consistency. Alternatively, relaxation effects should increase 
when the number of bodies decreases, and that was our sec- 
ond check. We took at random 10 per cent of the particles of 
each one of the models E5, increased their masses 10 times, 
and we run these new models self-consistently first for 150 
T cr to let them relax, and then for another 500 T cr to ana- 
lyze their stability. The corresponding changes are shown in 
Table [5] and, except for those of the central density which 
remain not significant, they are substantially larger than 
the equivalente ones of Table [3] Thus, we may conclude 
that even the small variations shown in the latter Table, 
are mainly due to relaxation effects of the N-body code. 

3.2 Regular, partially and fully chaotic orbits 

We randomly selected between 4, 500 and 5, 000 particles 
from each model and adopted their positions and velocities 
as the initial values to obtain the orbits and investigate their 
chaoticity. The potentials were fixed, keeping constant the 
coefficients of their expansions at their final values, and the 
integrations were carried out in fixed coordinate systems, in 
those cases where the rotation velocity was not significant 
(i.e., less than three times the mean square error), or in 
systems rotating with the corresponding velocity in those 
cases where it was significant. We proceeded in this way 
to be consistent with the models obtained, but even the 
significant velocities are so small that their effect is in all 
likelihood negligible (jMuzzidbOPl ). 

Since the potentials of our systems were fixed, all our 
orbits obey the energy integral (or the Jacobi integral in the 
case of systems with significant rotational velocity). There- 
fore, regular orbits have to obey two additional isolating 
integrals, but we can have two kinds of chaotic orbits: Par- 
tially chaotic orbits obey only one additional integral besides 
energy, and fully chaotic orbits have no isola ting integrals 
other than energy. We have shown before (see lMuzzio et ail 



120051 ; lAquilano et ai1l2007l ; iMuzzio et al]|2009l . for example) 
that partially and fully chaotic orbits have different distri- 
butions, so that it is important to separate them in orbital 
structure studies. A simple, albeit computing time demand- 
ing, way of classifying regular, partially and fully chaotic 
orbits is through the use of the six Lyapunov exponents. 
Since phase space volume is conserved, the exponents come 
in three pairs of the same absolute value and opposite sign. 
Due to energy conservation, one of those pairs is always zero 
in our case, and each additional isolating integral makes zero 
another pair. Thus, regular orbits have all their Lyapunov 
exponents equal to zero, partially chaotic orbits have one 
non-zero pair and fully chaotic orbits have two. 

The numerical equivalent of the Lyapunov exponents 
(which demand to integrate the orbit over an infinite 
time interval) are the finite time Lyapunov characteris- 
tic numbers (hereafter FT-LCNs) and, as in our previous 
works, we computed the m using the LIAMAG subroutine 
l|Udrv fc Pfennigerlll98Sh . kindly provided by D. Pfenniger, 
and that we adapted with the assistance of H.D. Navone 
to include the potential, accelerations and variational equa- 
tions corresponding to the method of Hernquist and Os- 
triker. The LIAMAG uses double precision arithmetic and a 
Runge - Kutta - Fehlberg integrator of order 7-8 with vari- 
able time steps. As in our previous works, we adopted the 
integration and normalization intervals as 10, 000 t.u. and 
lt.u., respectively, and we requested a precision of 10~ 15 for 
the step size regulation, which resulted in an energy conser- 
vation better than 2 x 10~ 12 at the end of integration in all 
cases. We will refer to the largest FT-LCN of a given orbit 
as L max and the second largest one as Lint, hereafter. 

Since the FT-LCNs are obtained from numerical inte- 
grations over a finite time interval (i.e., 10, 000 t.u. here), 
rather than the infinite one required to obtain Lyapunov 
exponents, they cannot reach zero value, but only a limiting 
minimum value, Lu m - Applying the definition of t h e Lya - 
punov exponents to equation (7) of ICincotta et all (|2003T ). 
one can estimate that value as Li im w InT/T, where T is 
the integration interval (Cincotta, private communication), 
that is, a limiting value of about 0.00092 (t.u.)^ 1 for our 
T — 10, 000 t.u. interval. This is an order of magnitude esti- 
mate only, however, and it is better to derive a more accu- 
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Table 4. Percentage variations over one Hubble time, for models E5 with constant coefficients. 



Mod Dens.(%) Xmom.in.(%) Ymom.in.{%) Zmom.in.(%) 

E5a +0.02 ± 0.74 -0.29 ± +0.09 -0.08 ± 0.10 +0.41 ± 0.07 

E5b +0.46 ± 0.80 -0.33 ± +0.06 -0.01 ± 0.10 +0.38 ± 0.13 

E5c -0.96 ± 0.44 -0.26 ± +0.05 -0.15 ± 0.12 +0.46 ± 0.10 



Table 5. Percentage variations over one Hubble time, for models E5 with ten times less bodies. 



Mod 


Dens.(%) 


Xmom.in.(%) 


Ymom.in.(%) 


Zmom.in.(%) 


E5a 


+1.18 ± 0.83 


-4.98 ± 0.50 


+5.70 ± 0.54 


+5.47 ± 0.53 


E5b 


+0.07 ± 1.11 


-4.75 ± 0.53 


+5.64 ± 0.45 


+6.38 ± 0.41 


E5c 


+0.68 ± 0.87 


-5.32 ± 0.46 


+5.21 ± 0.28 


+6.46 ± 0.76 



0.0025 



0.002 - 




Lmax 



Figure 4. Lmax versus L;„t for Model E5c. Only the region of 
low values, around the region corresponding to regular orbits is 
shown. 

rate value from the FT-LCNs themselves using, for example, 
plots of the low end of the L distribution (i.e., 

the region around that occupied by the representative points 
of the regular orbits) . As an example we show in Figure 0] 
the corresponding plot for model E5c. The representative 
points of the regular orbits are concentrated on the blob at 
the left and those of chaotic orbits extend towards the right 
(the sharp envelope of the identity line is merely due to the 
fact that, by definition, Li„t ^ Lmax). Clearly, any limit 
that attempts to separate regular from chaotic orbits can 
have a statistical value only, because some regular orbits 
may have FT-LCNs slightly larger than that limit, while 
some chaotic orbits may have somewhat lower FT-LCNs. 
With that caveat, using plots equivalent to that of Figure [4] 
for all our models we adopted Lu m — 0.0018 (t.u.)~ . 

Now, that Lnm corresponds to a Lyapunov time of 
556 t.u., which is equivalent to about 5 or 6 Hubble times 
for our models. Is it reasonable to use such a low Lu m to 
separate regular from chaotic orbits? Would it not be more 
sensible to adopt a value of 0.0100(t.u.) _1 which corresponds 
to a Lyapunov time of the same order of the Hubble time? 
Tho se questions have been considered in our previous work 
(see lAauilano et al ] |2007l :lM uzzio et al ] |2009l . for example) 
and we will repeat here the same analysis done there. 



First, we separated the orbits of each model into 
three groups: a) Those with Lmax < 0.0018 (t.u.)~ , 
i.e., those that are classified as regular for both 
choices of L; im (REGREG, hereafter); b) Those with 
0.0018 {t.u.)- 1 < L ma x < 0.0100 (t.u.) -1 , i-e., those 
that are classified as reg ular for L Hm = 0.0100 (t.u.)' 1 , but 
as chaotic for L iim = 0.0018 (t.u.)- 1 (REGCHAO, here- 
after); c) Those with 0.0100 (t. u. )" 1 ^ L max , ie, those that 
are classified as chaotic for both elections of Lu m (here- 
after CHAOCHAO). Then we considered, for each orbit, 11 
(x, y, z) orbital positions separated by intervals of 10 t.u., 
that is, over a total interval of 100 t.u., and, for each model 
and for each type of orbit, we computed the mean square 
value of each coordinate. Table [6] gives the square roots of 
the ratios of the y and z mean square values to the x mean 
square value. 

The results in Table IS] show that, as in our previous 
work, most of the axial ratios of the REGCHAO orbits are 
significantly different from those of the REGREG orbits and 
we may conclude that, despite their low FT-LCN values im- 
plying Lyapunov times longer than the Hubble time, orbits 
with 0.0018 (t.u.)" 1 ^ Lm ax < 0.0100 (t.u.)" 1 have a spatial 
distribution different from that of regular orbits. In other 
words, the time scale for the exponential divergence of or- 
bits (measured by the Lyapunov time) is not much rele- 
vant for the spatial distribution of those orbits. Since we are 
here interested in that distribution, the sensible approach is 
thus to adopt Lu m = 0.0018 (t.u.)" 1 . Therefore, we classify 
orbits as regular if L max < Lnm, as partially chaotic if 
Li n t < Lnm < Lmax and as fully chaotic if Lum < L int . 

Table [7] gives the percentages of regular, partially and 
fully chaotic orbits in our models. As in our previous works, 
the statistical errors have been estimated from the binomial 
distribution. The agreement among results of the a, b and 
c cases of each model is excellent and within the 3cr level 
in all cases, except for the regular and fully chaotic orbits 
of the E3c case where the differences with the E3a and E3b 
cases are of the order of 5<r. Clearly, all our models harbor 
the lowest fractions of regular orbits in triaxial systems we 
know about, l ower even than the 29.05 per cent found by 
iMuzzio et al.l |2009j) for their model E6c. Interestingly, the 
fraction of regular orbits does not decrease monotonically 
from type E2 through type E5, but falls first to rise again; 
a somewhat similar, but less p ronounced, trend was found 
in the E4 through E6 models of lAauilano et al.l (|2007r i. The 
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Table 6. Axial ratios of different classes of orbits for different 
choices of Lu m . 



Ratio Mod REGREG REGCHAO CHAOCHAO 



E2a 


1 


.030 


± 





015 





868 


± 





,013 





865 


± 





,015 


E2b 


1 


.023 


± 





014 





868 


± 





.013 





,875 


± 





,015 


E2c 


1 


.003 


± 





014 





847 


± 





,014 


0. 


856 


± 





,014 


E3a 





.816 


± 


0. 


022 


0. 


851 


± 





,018 





856 


± 





,018 


E3b 





.815 


± 





023 


0. 


824 


± 





,019 





829 


± 





,017 


E3c 





.773 


± 





025 





866 


± 





,018 


0. 


877 


± 





,016 


E4a 





.520 


± 


0. 


040 


0. 


745 


± 





,026 





874 


± 





019 


E4b 





.510 


± 





044 





788 


± 





,028 





841 


± 





,019 


E4c 





.568 


± 





039 





765 


± 





,023 


0. 


830 


± 





,019 


E5a 





.892 


± 


0. 


028 


0. 


877 


± 





,025 


0, 


908 


± 





018 


E5b 





.965 


± 





026 





833 


± 





,025 


0. 


899 


± 





,018 


E5c 





.934 


± 





029 





759 


± 





,025 





839 


± 





016 


E2a 





,910 


± 





016 





915 


± 





,013 





905 


± 





015 


E2b 





.937 


± 





015 





903 


± 





,013 





,874 


± 





,015 


E2c 





.879 


± 





015 





928 


± 





,014 





891 


± 





,014 


E3a 





.676 


± 


0. 


024 


0. 


840 


± 





,018 





832 


± 





,019 


E3b 





712 


± 





027 





812 


± 





,018 





771 


± 





,018 


E3c 





.668 


± 





028 


0. 


851 


± 





,018 


0. 


823 


± 





,017 


E4a 





.424 


± 


0. 


053 


0. 


719 


± 





,028 


0, 


760 


± 





,020 


E4b 





353 


± 





044 





660 


± 





,028 





767 


± 





,019 


E4c 





.401 


± 





049 





710 


± 





,026 


0. 


769 


± 





,020 


E5a 





.309 


± 


0. 


033 


0. 


501 


± 





,031 





759 


± 





021 


E5b 





361 


± 





040 





527 


± 





,035 





699 


± 





,020 


E5c 





.319 


± 





043 





473 


± 





,032 





699 


± 





,018 



decreasing fractions of partially chaotic orbits when going 
from type E2 through E5 show the same trend found in 
those two previous works of ours. 

Two further tests were performed to check the accuracy 
of our results. On the one hand, it is interesting to check the 
influence that the use of a particular N-body snapshot as a 
model might have. On the other hand, since we are using 
only about 0.5 per cent of the orbits in each model to inves- 
tigate chaoticity, it might be worthwile to see if the use of a 
different sample yields different results. We adopted the E4b 
model for both tests. For the first one, we took a new model 
from the snapshot evolved 200T cr more than the one whose 
results are presented in Table [3 and we obtained fractions 
of 14.42 ± 0.52%, 11.95 ± 0.48% and 73.62 ± 0.65%, respec- 
tively for the regular, partially and fully chaotic orbits. The 
first two fractions coincide with those of the original model 
at the 3cr level, but the last one differs by 3.4<r. Not only are 
these differences small, but they are probably at least par- 
tially due to the axial changes of the models as they evolve 
described in Subsection 3.1. The 200T cr evolution changes 
the triaxiality of the E4b model in the direction of that of 
the E5 models by about 1/10 of the triaxiality difference 
between the E4 and E5 models. Since the difference be- 
tween the fractions of fully chaotic orbits of those models 
is about 8% we might crudely estimate that about 0.8% of 
the change from the original model to the one evolved 200T cr 
might be attributed to the axial changes, in which case the 



Table 7. Percentages of regular and chaotic orbits in triaxial 
systems. 



Mod Regular(%) Part. Chaotic(%) Fully Chaotic(%) 

E2a 
E2b 
E2c 

E3a 
E3b 
E3c 

E4a 
E4b 
E5c 

E5a 
E5b 
E5c 



remaining difference is of 2.5<r only. Nevertheless, the agree- 
ment between the original and the 200T cr evolution results 
is somewhat poorer than those among statistically different 
realizations of the same model, except for the already men- 
tioned differences for models E3. The second test, instead, 
yielded no significant differences. The new sample of orbits 
from E4b model had fractions of 12.50±0.49%, 9.74±0.44% 
and 77.76 ±0.62%, respectively of regular, partially and full 
chaotic orbits, all of them well within the 3cr level of differ- 
ences from the original results. 

The values of the FT-LCNs are larger than those we 
fo und for non-cuspy models , in agreement with the result s 
of iKandrup fc Siderisl (|2002h and lKandrup fc Siopisl l|2003h . 
Figure [S] shows, for our model E5c, the plot of L max ver- 
sus the reduced energy (i.e., the orbital energy, E, divided 
by the potential energy at the galactic centre, W ) , which 
can be compared with Fig. 3 of iMuzzio et al. | (|2005h . In the 
present case the largest FT-LCN s are close to 1.0 (t. u.)~ , 
while for the non-cuspy model of IMuzzio et al.l (|2005h they 
only reached about 0.5 (t.u.) -1 . Besides, in the non-cuspy 
system there were no chaotic orbits for E/W values close 
to 1.0 (i.e., orbits very close to the centre of the system), 
because the potential near the centre was essentially that of 
a three dimensional harmonic oscillator. But, here, the pres- 
ence of the cusp strongly limits the central region where the 
potential can be approximated with an integrable potential 
and, therefore, one can find chaotic orbits near the centre of 
the system. 

For each model we took eleven orbital positions at in- 
tervals of 10 t.u., i.e., over a total interval of 100 t.u., to 
obtain the mean square values of each coordinate separately 
for each type of orbit. Table [8] gives the y/x and z/x axial 
ratios, computed from the square roots of those quadratic 
mean values. All the results for the a, b and c realizations 
of a same model show very good agreement at the 3<r level, 
except for the y/x ratio difference between models E4b and 
E4c for the fully chaotic orbits, which is somewhat larger. 
The differences between the distributions of the partially 
and fully chaotic orbits are significant at the 3er level only 
for z/x ratio of the E5 models and the b case of models 
E4. This result a grees with those obtained for the models 
we studied before (|Aquilano et al.ll2007l ; IMuzzio et~al]|2009h 



22 


.48 


± 





,59 


21 


.35 


± 





,58 


22 


.24 


± 





,59 


11 


.21 


± 





50 


11 


.04 


± 





50 


10 


.63 


± 





11 


13 


,05 


± 





,50 


12 


72 


± 





50 


12 


.67 


± 





19 


21 


.71 


± 





61 


23 


,06 


± 





.63 


20 


.92 


± 





.60 



15 


.30 


± 


0.51 


15 


,52 


± 


0.51 


13 


,64 


± 


0.49 


13 


19 


± 


0.48 


13 


,57 


± 


0.49 


12 


,90 


± 


0.48 


10 


,85 


± 


0.46 


10 


,58 


± 


0.46 


10 


11 


± 


0.45 



9.69 ± 0.44 
9.38 ± 0.43 
8.86 ± 0.42 



62 


22 


± 


0. 


69 


63 


.13 


± 


0. 


69 


61 


12 


± 





68 


72 


60 


± 


0. 


64 


72 


.39 


± 


0. 


64 


76 


,17 


± 


0. 


61 


76 


10 


± 


0. 


63 


76 


70 


± 


0. 


63 


77 


22 


± 





62 


68 


60 


± 


0. 


69 


67 


.96 


± 


0. 


69 


70 


22 


± 





68 
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4 CONCLUSIONS 



E 0.5 




Figure 5. The maximum FT-LCNs (L max ) versus reduced or- 
bital energy plot for Model E5c 



Table 8. Axial ratios of different kinds of orbits. 



Ratio 


Mod 


Regular 


Part. Chaotic 


Fully Chaotic 


y/x 


E2a 


1.030 ± 


0.015 


0.903 ± 0.016 


0.828 ± 0.014 




E2b 


1.023 ± 


0.014 


0.899 ± 0.016 


0.841 ± 0.014 




E2c 


1.003 ± 


0.014 


0.891 ± 0.016 


0.813 ± 0.014 




E3a 


0.816 ± 


0.022 


0.878 ± 0.020 


0.800 ± 0.024 




E3b 


0.815 ± 


0.023 


0.863 ± 0.021 


0.765 ± 0.027 




E3c 


0.773 ± 


0.025 


0.915 ± 0.020 


0.799 ± 0.025 




E4a 


0.520 ± 


0.040 


0.804 ± 0.032 


0.712 ± 0.032 




E4b 


0.510 ± 


0.044 


0.781 ± 0.035 


0.814 ± 0.034 




E4c 


0.568 ± 


0.039 


0.894 ± 0.030 


0.652 ± 0.029 




E5a 


0.892 ± 


0.028 


0.869 ± 0.029 


0.903 ± 0.018 




E5b 


0.965 ± 


0.026 


0.868 ± 0.030 


0.871 ± 0.017 




E5c 


0.934 ± 


0.029 


0.790 ± 0.029 


0.813 ± 0.016 


z/x 


E2a 


0.910 ± 


0.016 


0.905 ± 0.016 


0.919 ± 0.013 




E2b 


0.937 ± 


0.015 


0.870 ± 0.016 


0.917 ± 0.013 




E2c 


0.879 ± 


0.015 


0.930 ± 0.016 


0.903 ± 0.014 




E3a 


0.676 ± 


0.024 


0.825 ± 0.021 


0.867 ± 0.024 




E3b 


0.712 ± 


0.027 


0.791 ± 0.021 


0.831 ± 0.026 




E3c 


0.668 ± 


0.028 


0.845 ± 0.022 


0.850 ± 0.025 




E4a 


0.424 ± 


0.053 


0.648 ± 0.037 


0.798 ± 0.034 




E4b 


0.359 ± 


0.045 


0.573 ± 0.037 


0.789 ± 0.033 




E4c 


0.401 ± 


0.049 


0.699 ± 0.036 


0.730 ± 0.031 




E5a 


0.309 ± 


0.033 


0.374 ± 0.033 


0.733 ± 0.020 




E5b 


0.361 ± 


0.040 


0.484 ± 0.045 


0.675 ± 0.019 




E5c 


0.319 ± 


0.043 


0.398 ± 0.038 


0.671 ± 0.018 



in the sense that the differences between the distributions 
of partially and fully chaotic orbits, on the one hand, in- 
creased for more flattened systems and, on the other hand, 
diminished for cuspy systems. 



The method of Hcrnquist & Ostrikcr ( 1992) is clearly better 
than that of Aguilar ( Aguilar fc Merritt] 1990t ) to build cuspy 
triaxial stellar models, since the former needs no softening 
and uses a radial expansion which is particularly adequate 
to fit N-body distributions with 7 ~ 1. As a side benefit, 
the coefficients of both the radial and angular expansions 
are provided by the method itself and there is no need of 
additional fit tings to obtain them. The on ly caveat is that, 
as shown by iKalapotharakos et alj (|200ST > . the method of 
Hernquist and Ostriker is adequate for cases with 7 ~ 1 only 
and, for other slopes of the cusp, other related but different 
methods should be preferred. 

The models obtained using the same parameters and 
changing only the seed number of the random number gen- 
erator show that our results are very robust, indeed. Those 
models exhibit essentially the same global properties (ex- 
cept for the angular velocity of figure rotation) and, with 
very few exceptions, the same fractions and distributions 
of regular, partially and fully chaotic orbits. The different 
angular velocities we obtained here for statistically equiva- 
lent models further complicates the little we know about the 
phenomenon of figure rotation (i.e., with zero angular mo- 
mentum) in this kind of models. On the one hand, the real- 
ity of the rotation has b een checked us ing both the A guilar 
llAguilar fc Merritt|[T990l ) and Aarseth (|Aarsetlj|20"03l ) codes 
bv lMuzzid I (|2006h and here with the Hernquist code and, on 
the other hand, such triaxial systems with figure rotation are 
just the stellar equivalent of the Riemann ellipsoids of fluid 
dynamics, so that there is nothing misterious about their ex- 
istence. Nevertheless, except for the increase of the angular 
velocity when one goes to flatter systems, the rotation does 
not seem to correlate with other properties of the systems 
and, moreover, now it turns out that statistically equivalent 
systems have different velocities. Even though, as already 
indicated, the angular velocities are so low that they have 
very little effect on orbital chaoticity, this phenomenon is 
interesting and warrants further investigation. 

The high stability of our models over time scales of the 
order of a Hubble time is very well stablished. Not only are 
the detected variations very small, but they are most likely 
due to relaxation effects of the N-body code. Besides, it 
should be stressed that we have checked the stability running 
the models self-consistently, while che cks on the stab i lity o f 
models resulting from the method of ISchwarzschildl (Il979h 
are usually done on fixed potentials ()Schwarzschildl fl993l ) . 
As we have shown, when the potential is fixed, the changes 
in our models are an order of magnitude smaller than those 
found self-consistently. 

Al though, considering the re sults of Kandrup fc Siderisl 

l|2002t l. iKandrup fc Siopisf (|2003h and lMuzzio et al.l l|2009h . 
high fractions of chaotic orbits could be expected in our 
cuspy triaxial models, the fractions found are extremely 
high, indeed, exceeding 85 per cent for our E2 and E3 mod- 
els. It is important to emphasize that those high fractions of 
chaotic orbits are not merely due to the low Lu m adopted: 
Had we adopted the 0.0T0(t.u.) -1 limit, those models would 
still have had between 70 an 75 per cent of chaotic orbits. 

The main conclusion of our investigation is that it 
is perfectly possible to build self-consistent stellar mod- 
els of cuspy triaxial elliptical galaxies that are very sta- 
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ble despite containing high percentages of chaotic orbits. 
Thi s result confirms and extends previous work by oth- 
ers llVoglis et al. II2OO2I: iKalapotharakos fc Voglis! 12009 ) and 
ours (|Aquilano et al.ll2007l; iMuzzio et alj|2009h . It also sup- 
ports the suggestion by IMuzzio et al, ( 20051 ) in the sense 
that the difficulties to obtain such models with the method 
of Schwarzschild should be attributed to the method itself 
and not to physical causes. It is true that, while regular or- 
bits always occupy the same region of space, chaotic orbits 
may spend a long interval occupying certain region, only 
to switch to a different one later on, with sticky orbits be- 
ing an extreme example. But, although such behaviour con- 
spires against obtaining stable models with high fractions 
of chaotic orbits with the method of Schwarzschild, it does 
not necessarily prevents the existence of those models. All 
one needs is that, as some chaotic orbits abandon a region 
of space to explore a different one, orbits in the latter re- 
gion move on to replace those in the former. In other words, 
rather than having a static equilibrium, with each orbit cov- 
ering always the same zone, we have a dynamic equilibrium 
where switches from one zone to another are present but 
balance each other on average. Our models show that such 
state of affairs is perfectly possible. 
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